# This script runs all the analyses reported in figures and tables 
# in the main text and in the appendix

# 2.4. Measuring Support for Putin ----------
# Figure 1 ----------------
# the distribution of presidential approval across three surveys
pa_means <- map2(list(main_resp |>
                        filter(session_id %in% main_stories$session_id),
                      levada_resp,
                      omi_resp),
                 list("Main\nstudy",
                      "National",
                      "Online\npanel"),
                 ~.x |>
                   drop_na(pres_approval_cat) |>
                   rename(approval = pres_approval_cat) |>
                   group_by(approval) |>
                   summarise(share = n()) |>
                   ungroup() |>
                   mutate(share = round(share/sum(share)*100, 1),
                          sample = .y)) |>
  list_rbind() |>
  mutate(approval = str_replace(approval, " ", "\n"),
         approval = factor(approval, 
                           levels = c("Certainly\napprove",
                                      "Somewhat\napprove",
                                      "Somewhat\ndisapprove",
                                      "Certainly\ndisapprove")),
         sample = factor(
           sample, 
           levels = c("Main\nstudy",
                      "Online\npanel", 
                      "National")))

# save estimates for the paper
write_csv(pa_means, "Estimates/Figure_1_approval_dist_3_surveys.csv")

# Figure 1
ggplot(pa_means, aes(sample, share, color = sample, fill = sample)) +
  geom_bar(stat = "identity", width = 0.9, 
           position = position_dodge(width = 0.1)) +
  geom_text(aes(label = share), 
            position = position_stack(vjust = 0.5),
            color = "white") +
  scale_color_manual(values = c("#1b9e77", "#d95f02", "#7570b3")) +
  scale_fill_manual(values = c("#1b9e77", "#d95f02", "#7570b3")) +
  labs(x = "", y = "", 
       title = "Distribution of Putin's approval in the surveys, %") +
  theme(legend.position = "none") +
  facet_wrap(~approval, ncol = 2) +
  theme_classic() +
  theme(legend.position = "none")

ggsave("Figures/Figure_1.png")

# 3.1. Supporters Are Receptive to Identity-Consistent Stories... -----------
# Regressions for Figures 2 and B2 --------
# regression for Figure 2 (and left panel of Figure B2), 
# stories labeled "Pro-Russia stories, mean" and "Critical stories, mean"
stories_perc <- 
  reg_run(c("pres_approval_dummy",
            "story_direction", "story_current",
            "female", "age", "education",
            "took_quiz_2", 
            "story_order"),
          dep_var = "story_real",
          interaction_terms = c("pres_approval_dummy",
                                "story_direction"),
          story_set = "all",
          first_three_stories = T,
          contrast_var = "pres_approval_dummy",
          contrast_filter = "pres_approval_dummy0 - pres_approval_dummy1")

# regression for Figure 2, remaining story categories ("Praising Putin", etc.)
stories_perc_detailed <- 
  reg_run(c("pres_approval_dummy",
            "story_subdirection", "story_current",
            "female", "age", "education",
            "took_quiz_2", 
            "story_order"),
          dep_var = "story_real",
          interaction_terms = c("pres_approval_dummy",
                                "story_subdirection"),
          story_set = "all",
          first_three_stories = T,
          contrast_var = "pres_approval_dummy",
          contrast_filter = "pres_approval_dummy0 - pres_approval_dummy1")

# regression for the central panel of Figure B2 
# ("All supporters vs critics, no sources")
stories_perc_nos <- 
  reg_run(c("pres_approval_dummy",
            "story_direction", "story_current",
            "female", "age", "education",
            "took_quiz_2", 
            "story_order"),
          dep_var = "story_real",
          interaction_terms = c("pres_approval_dummy",
                                "story_direction"),
          story_set = "all",
          first_three_stories = T,
          no_source_stories = T,
          contrast_var = "pres_approval_dummy",
          contrast_filter = "pres_approval_dummy0 - pres_approval_dummy1")

# regression for the right panel of Figure B2 
# ("Strong vs moderate supporters")
stories_perc_supporters <- 
  reg_run(c("pres_approval_cat",
            "story_direction", "story_current",
            "female", "age", "education",
            "took_quiz_2", 
            "story_order"),
          dep_var = "story_real",
          interaction_terms = c("pres_approval_cat",
                                "story_direction"),
          story_set = "all",
          first_three_stories = T,
          contrast_var = "pres_approval_cat",
          contrast_filter = c("Somewhat approve - Certainly approve"))

# combine contrasts in one data frame for plotting figures
stories_perc_all <- bind_rows(stories_perc$effect |>
                                mutate(contrast = 
                                         "All supporters vs critics",
                                       model = "Average"),
                              stories_perc_nos$effect |>
                                mutate(contrast = 
                                         "All supporters vs critics,\nno sources",
                                       model = "Average"),
                              stories_perc_supporters$effect |>
                                mutate(contrast = 
                                         "Strong vs moderate\nsupporters",
                                       model = "Average"),
                              stories_perc_detailed$effect |>
                                mutate(contrast = 
                                         "All supporters vs critics",
                                       model = "Detailed") |>
                                rename(story_direction = story_subdirection)) |>
  select(contrast, story_direction, estimate, conf.low, conf.high, model) |>
  rename(Stories = story_direction, diff = estimate) |>
  # transform variables/additional variables for plots
  mutate(across(c(diff, conf.low, conf.high), ~100*.x),
         `Story content` = case_when(
           Stories %in% c("Praising Ukraine", 
                          "Criticizing Russia",
                          "Criticizing Putin", 
                          "Critical") ~ "Critical",
           Stories %in% c("Praising Putin", 
                          "Praising Russia",
                          "Criticizing West",
                          "Criticizing Ukraine",
                          "Pro-Russia") ~ "Pro-Russia",
           TRUE ~ "Other"
         ),
         Stories = if_else(Stories %in% c("Critical", "Pro-Russia"),
                           paste0(Stories, " stories, mean"),
                           Stories)) |>
  filter(`Story content` != "Other")

# save estimates for the paper (Figures 2 and B2)
write_csv(stories_perc_all, "Estimates/Figure_2_B2_stories_perc_by_approval_adj.csv")

# Figure 2 ------
ggplot(stories_perc_all |> 
         filter(contrast == "All supporters vs critics") |>
         mutate(Stories = factor(Stories, 
                                 levels = c(
                                   "Praising Ukraine", "Criticizing Russia",
                                   "Criticizing Putin", "Critical stories, mean", 
                                   "Criticizing West", "Criticizing Ukraine",  
                                   "Praising Russia", "Praising Putin", 
                                   "Pro-Russia stories, mean"))), 
       aes(Stories, diff, linetype = `Story content`, shape = `Story content`, 
           color = model, alpha = model)) +
  geom_point(size = 3,
             position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(ymin = conf.low,
                    ymax = conf.high),
                width = 0,
                position = position_dodge(width = 0.4),
                size = 1.1) +
  geom_text(aes(label = round(diff, 1)), vjust = -0.5, size = 4, 
            show.legend = F) +
  scale_color_manual(values = c("black", "darkgray"), guide = "none")  +
  scale_alpha_manual(values = c(0.6, 1), 
                     guide = "none") +
  geom_hline(yintercept = 0, linetype = "dashed", col = "gray") +
  labs(x = "", 
       y = "Difference in the probability of saying\nstory is true, percentage points, 95% CI",
       title = "Belief in pro-Russia and critical stories",
       subtitle = "Difference between Putin supporters and opposition-minded respondents") +
  ylim(c(-20, 20)) +
  coord_flip() +
  theme_classic() +
  theme(legend.position = "bottom")
ggsave("Figures/Figure_2.png")

# Figure B2 --------
ggplot(stories_perc_all |> 
         filter(model == "Average"), 
       aes(Stories, diff)) +
  geom_point(size = 3,
             position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(ymin = conf.low,
                    ymax = conf.high),
                width = 0,
                position = position_dodge(width = 0.4),
                size = 1) +
  geom_text(aes(label = round(diff, 1)), vjust = -0.5, size = 4.5, 
            show.legend = F) +
  geom_hline(yintercept = 0, linetype = "dashed", col = "gray") +
  labs(x = "", y = "Difference, percentage points",
       title = "") +
  coord_flip() +
  ylim(c(-20, 20)) +
  facet_wrap(~contrast, ncol = 3) +
  theme_classic() +
  theme(legend.position = "bottom")

ggsave("Figures/Figure_B2.png")

# Figure B3 -------
# regressions
# main study
# run a model to calculate the contrast between supporters and critics
# for every single story
stories_perc_selected_main <- 
  reg_run(c("pres_approval_dummy",
            "story_label",
            "female", "age", "education",
            "story_order"),
          dep_var = "story_real",
          interaction_terms = c("pres_approval_dummy",
                                "story_label"),
          story_set = "all",
          first_three_stories = T,
          story_dummies = F,
          contrast_var = "pres_approval_dummy",
          contrast_filter = "pres_approval_dummy0 - pres_approval_dummy1")

# select only stories that are also included in the OMI and Levada surveys
stories_perc_selected_main_contrast <- stories_perc_selected_main$effect |>
  filter(as.character(story_label) %in% c("103", "107", "110", "111")) |>
  mutate(story_label = as.character(story_label),
         story_label = case_when(
           story_label == "103" ~ "Growth in Ukraine",
           story_label == "107" ~ "Trump and Putin",
           story_label == "110" ~ "COVID origins",
           story_label == "111" ~ "Nuclear waste"
         )) 

# Levada study
stories_perc_selected_levada <- lm_robust(story_response_scaled ~ pres_approval_dummy  +
                                            story_label + female + age + education +  
                                            pres_approval_dummy*story_label +
                                            story_source_tv1,
                                          data = levada_stories,
                                          weights = weight,
                                          clusters = respondent_id,
                                          se_type = "CR0")
# contrast for every story
stories_perc_selected_levada_emm <- emmeans(stories_perc_selected_levada, 
                                            ~pres_approval_dummy*story_label)
stories_perc_selected_levada_contrast <- summary(pairs(
  stories_perc_selected_levada_emm, simple = "pres_approval_dummy")) |>
  # select the story that was also in the main study
  filter(story_label %in% c("Ukrainian economy vs Russian")) |>
  # reverse the sign so that the difference is "supporters" - "critics"
  # as in the analysis of the main study
  mutate(estimate = -estimate,
         conf.low = estimate - 1.96*SE,
         conf.high = estimate + 1.96*SE,
         story_label = as.character(story_label),
         story_label = "Growth in Ukraine")

# OMI study
stories_perc_selected_omi <- lm_robust(story_real ~ pres_approval_dummy  +
                                         story_label + female + age + education +  
                                         pres_approval_dummy*story_label +
                                         story_source,
                                       data = omi_stories |>
                                         mutate(age = age_ordinal,
                                                education = education_higher),
                                       clusters = ResponseId,
                                       se_type = "CR0")
# contrast for every story
stories_perc_selected_omi_emm <- emmeans(stories_perc_selected_omi, 
                                         ~pres_approval_dummy*story_label)
stories_perc_selected_omi_contrast <- summary(pairs(stories_perc_selected_omi_emm, 
                                                    simple = "pres_approval_dummy")) |>
  # select the story that was also in the main study
  filter(story_label != "mafia") |>
  # reverse the sign so that the difference is "supporters" - "critics"
  # as in the analysis of the main study
  mutate(estimate = -estimate,
         conf.low = estimate - 1.96*SE,
         conf.high = estimate + 1.96*SE,
         story_label = case_when(
           story_label == "trump" ~ "Trump and Putin",
           story_label == "uranium" ~ "Nuclear waste",
           story_label == "virus" ~ "COVID origins",
         ))

# combine all contrasts for the plot
stories_perc_selected_all <- bind_rows(stories_perc_selected_main_contrast |>
                                         mutate(Survey = "Main study"),
                                       stories_perc_selected_omi_contrast |>
                                         mutate(Survey = "Online panel"),
                                       stories_perc_selected_levada_contrast |>
                                         mutate(Survey = "National")) |>
  mutate(story_direction = case_when(
    story_label %in% c("Growth in Ukraine",
                       "Nuclear waste") ~ 
      "Critical stories",
    TRUE ~ "Pro-Russia stories"
  ),
  across(c(estimate, conf.low, conf.high), ~100*.x))

# save estimates for the paper
write_csv(stories_perc_selected_all,
          "Estimates/Figure_B3_stories_perc_by_approval_3_surveys.csv")

# Figure B3
ggplot(stories_perc_selected_all, 
       aes(story_label, estimate, linetype = Survey, shape = Survey)) +
  geom_point(size = 2,
             position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(ymin = conf.low,
                    ymax = conf.high),
                width = 0,
                position = position_dodge(width = 0.4),
                size = 0.8) +
  geom_hline(yintercept = 0, linetype = "dashed", col = "gray") +
  labs(x = "", y = "Difference between\nsupporters and critics",
       title = "Supporters and critics disagree about news stories in all samples") +
  ylim(c(-30, 27)) +
  facet_wrap(~story_direction, ncol = 2, scales = "free_x") +
  theme_classic() +
  theme(legend.position = "bottom")

ggsave("Figures/Figure_B3.png")

# 3.2 Supporters Find Propaganda Outlets More Credible... ------------
# Regressions for Figure 3 and Tables B4 and B8 -----
# model for the main treatment effect: black circled estimates in Figure 3
# and "Model 1" in Table B4
main_treatment_effect <- 
  reg_run(c("state_controlled", 
            "pres_approval_cat", "story_order"), 
          interaction_terms = c("state_controlled", 
                                "pres_approval_cat"),
          story_dummies = T)

# model for the treatment effect conditional on story direction:
# the remaining estimates in Figure 3 and "Model 1" in Table B8
main_treatment_effect_by_direction <- 
  reg_run(c("state_controlled", 
            "pres_approval_cat", "story_order", "story_direction"), 
          interaction_terms = c("state_controlled", 
                                "pres_approval_cat", 
                                "story_direction"),
          story_dummies = F)

# main treatment effect, adjusted for sociodemographic covariates
# "Model 2" in Table B4
main_treatment_effect_adj <- 
  reg_run(c("state_controlled", "age", "female", "education",
            "pres_approval_cat", "story_order"), 
          interaction_terms = c("state_controlled", 
                                "pres_approval_cat"),
          story_dummies = T,
          quantities_calculate = "model only")

# main treatment effect, logit model
# "Model 3" in Table B4
main_treatment_effect_logit <- 
  reg_run(c("state_controlled",
            "pres_approval_cat", "story_order"), 
          interaction_terms = c("state_controlled", 
                                "pres_approval_cat"),
          story_dummies = T,
          logit_model = T)

# effect estimates for Figure 3
# combine the effects conditional on presidential approval
# and conditional on approval and story content
main_treatment_effect_estimates <-
  bind_rows(
    main_treatment_effect$effect |>
      mutate(story_direction = "All stories"),
    main_treatment_effect_by_direction$effect) |>
  mutate(story_direction = factor(story_direction, 
                                  levels = c("Critical",
                                             "Pro-Russia", 
                                             "Neutral",
                                             "All stories")),
         pres_approval_cat = case_when(
           pres_approval_cat == "Certainly disapprove" ~ "Strong\ncritic",
           pres_approval_cat == "Somewhat disapprove" ~ "Moderate\ncritic",
           pres_approval_cat == "Somewhat approve" ~ "Moderate\nsupporter",
           pres_approval_cat == "Certainly approve" ~ "Strong\nsupporter"
         ),
         pres_approval_cat = factor(pres_approval_cat, 
                                    levels = c("Strong\ncritic",
                                               "Moderate\ncritic",
                                               "Moderate\nsupporter",
                                               "Strong\nsupporter")),
         across(c(estimate, conf.low, conf.high), ~100*.x))

# save estimates for the paper
write_csv(main_treatment_effect_estimates, 
          "Estimates/Figure_3_main_treatment_effect_estimates.csv")

# Figure 3 -----
ggplot(main_treatment_effect_estimates, 
       aes(pres_approval_cat, estimate,
           shape = story_direction, linetype = story_direction, 
           color = story_direction)) +
  geom_point(size = 3,
             position = position_dodge(width = 0.6)) +
  geom_errorbar(aes(ymin = conf.low,
                    ymax = conf.high),
                width = 0,
                position = position_dodge(width = 0.6),
                size = 1.1) +
  scale_color_manual(values = c("gray25", "gray50", "gray75", "black")) +
  scale_linetype_manual(values = c("dashed", "dotted", "dotted", "solid")) +
  scale_shape_manual(values = c(18, 15, 17, 19)) +
  geom_hline(yintercept = 0, linetype = "dashed", col = "gray") +
  labs(x = "", y = "Effect of switching the source from independent to state media",
       title = "Supporters find state media more credible",
       subtitle = "(except when critical stories are considered)",
       linetype = "Story group", shape = "Story group", color = "Story group") +
  ylim(c(-15, 15)) +
  coord_flip() +
  theme_classic() +
  theme(legend.position = "bottom")
ggsave("Figures/Figure_3.png")


# APPENDIX Table B4 ----------
texreg(list("Model 1" = main_treatment_effect$model,
            "Model 2" = main_treatment_effect_adj$model,
            "Model 3" = main_treatment_effect_logit$model),
       digits = 3,
       omit.coef = "survey_day|story_label",
       custom.coef.names = c("Intercept", "Source: Critical", 
                             "Source: State-controlled",
                             "Source: RBC", 
                             "Somewhat disapprove", "Somewhat approve",
                             "Certainly approve", "Story order", 
                             "Source: Critical*Somewhat disapprove",
                             "Source: State-controlled*Somewhat disapprove",
                             "Source: RBC*Somewhat disapprove",
                             "Source: Critical*Somewhat approve",
                             "Source: State-controlled*Somewhat approve",
                             "Source: RBC*Somewhat approve",
                             "Source: Critical*Certainly approve",
                             "Source: State-controlled*Certainly approve",
                             "Source: RBC*Certainly approve",
                             "Age",
                             "Female",
                             "Higher education"),
       label = "tabonlinestatemediaeffect",
       caption = paste("Treatment effect in the main study"),
       custom.note = paste("\\item %stars. Estimates from regression models",
                           "(OLS in Models 1 and 2, logit in Model 3)",
                           "with news story evaluations as dependent variables.",
                           "The reference category in presidential approval",
                           "is 'Certainly disapprove.'",
                           "The reference category in source treatments",
                           "is 'No source.'",
                           "Data from the social media sample.",
                           "Story and day fixed effects included.",
                           "Heteroskedasticity-robust standard errors",
                           "in parentheses (clustered on respondent)."),
       #       single.row = T,
       fontsize = "scriptsize",
       float.pos = "H",
       caption.above = TRUE,
       threeparttable = TRUE,
       use.packages = FALSE,
       file = "Tables/Table_B4_main_treatment_effect.tex",
       include.ci = FALSE)

# APPENDIX Figure B4 and Table B7  --------
# Other Measures of Pro-Regime Orientations
# a version of the model in Figure 3 where presidential approval is replaced 
# with an alternative measure based on beliefs about EU sanctions/Ukraine economy
main_treatment_effect_eu_ukr <- 
  reg_run(c("state_controlled",
            "EU_Ukr_beliefs", "story_order"), 
          interaction_terms = c("state_controlled", 
                                "EU_Ukr_beliefs"),
          story_dummies = T)

# a version of the model in Figure 3 where presidential approval is replaced 
# with an alternative measure, the dummy indicating pride in Crimea annexation
main_treatment_effect_pride_crimea <- 
  reg_run(c("state_controlled",
            "pride_history_crimea_cat", "story_order"), 
          interaction_terms = c("state_controlled", 
                                "pride_history_crimea_cat"),
          story_dummies = T)

# combine estimates for the paper/Figure B4
main_treatment_effect_eu_ukr_est <- 
  main_treatment_effect_eu_ukr$effect |>
  mutate(across(c(estimate, conf.low, conf.high), ~100*.x))
main_treatment_effect_pride_crimea_est <- 
  main_treatment_effect_pride_crimea$effect |>
  mutate(across(c(estimate, conf.low, conf.high), ~100*.x),
         pride_history_crimea_cat = 
           if_else(str_detect(pride_history_crimea_cat, "Not proud"), 
                   "Not proud", "Proud of Crimea"))

# save estimates for the paper
write_csv(main_treatment_effect_eu_ukr_est, 
          "Estimates/Figure_B4_main_treatment_effect_eu_ukr.csv")
write_csv(main_treatment_effect_pride_crimea_est, 
          "Estimates/Figure_B4_main_treatment_effect_pride_crimea.csv")

# Figure B4
# left panel
b4left <- ggplot(main_treatment_effect_pride_crimea_est, 
                 aes(pride_history_crimea_cat, estimate)) +
  geom_point(position = position_dodge(width = 0.4), size = 2) +
  geom_errorbar(aes(ymin = conf.low,
                    ymax = conf.high),
                width = 0,
                position = position_dodge(width = 0.4),
                size = 0.8) +
  geom_text(aes(label = round(estimate, 1)), vjust = -0.5, size = 3, 
            show.legend = F) +
  geom_hline(yintercept = 0, linetype = "dashed", col = "gray") +
  labs(x = "", y = "Effect of state media",
       title = "Given feelings toward Crimea") +
  ylim(c(-17, 17)) +
  theme_classic() +
  coord_flip()

# right panel
b4right <- ggplot(main_treatment_effect_eu_ukr_est, 
                  aes(EU_Ukr_beliefs, estimate)) +
  geom_point(position = position_dodge(width = 0.4), size = 2) +
  geom_errorbar(aes(ymin = conf.low,
                    ymax = conf.high),
                width = 0,
                position = position_dodge(width = 0.4),
                size = 0.8) +
  geom_text(aes(label = round(estimate, 1)), vjust = -0.5, size = 3, 
            show.legend = F) +
  geom_hline(yintercept = 0, linetype = "dashed", col = "gray") +
  labs(x = "", y = "Effect of state media",
       title = "Given beliefs about EU/Ukraine") +
  ylim(c(-17, 17)) +
  theme_classic() +
  coord_flip()

b4left + b4right + plot_annotation(title = "Effect of state vs critical media")
ggsave("Figures/Figure_B4.png")

# Table B7
texreg(list("Model 1" = main_treatment_effect_pride_crimea$model,
            "Model 2" = main_treatment_effect_eu_ukr$model),
       digits = 3,
       omit.coef = "survey_day|story_label",
       custom.coef.names = c("Intercept", "Source: Critical", "Source: State-controlled",
                             "Source: RBC", "Proud of Crimea",
                             "Story order", 
                             "Source: Critical*Proud of Crimea",
                             "Source: State-controlled*Proud of Crimea",
                             "Source: RBC*Proud of Crimea",
                             "EU-Ukraine beliefs: In-between",
                             "EU-Ukraine beliefs: Pro-regime",
                             "Critical*EU-Ukraine In-between",
                             "State-controlled*EU-Ukraine In-between",
                             "Source: RBC*EU-Ukraine In-between",
                             "Source: Critical*EU-Ukraine Pro-regime",
                             "Source: State-controlled*EU-Ukraine Pro-regime",
                             "Source: RBC*EU-Ukraine Pro-regime"),
       label = "tabonlinestatemediaeffectaltapproval",
       caption = paste("Treatment effect in the main study",
                       "(alternative measures of pro-regime attitudes)"),
       custom.note = paste("\\item %stars. Estimates from linear regressions",
                           "with news story evaluations as dependent variables.",
                           "Data from the social media sample.",
                           "Approval measures: pride in Crimea annexation",
                           "(Model 1), beliefs about EU and Ukraine",
                           "(Model 2); see text for details",
                           "Story and day fixed effects included.",
                           "Heteroskedasticity-robust standard errors",
                           "in parentheses (clustered on respondent)."),
       #       single.row = T,
       float.pos = "H",
       fontsize = "scriptsize",
       caption.above = TRUE,
       threeparttable = TRUE,
       use.packages = FALSE,
       file = "Tables/Table_B7_main_treatment_effect_alt_def.tex",
       include.ci = F)

# APPENDIX Figure B5 and Table B5  ---------------- 
# experimental results by individual sources (news outlets)
# model where presidential approval is interacted with the assigned news source
# (instead of the state-controlled dummy that was used for Figure 3)
main_treatment_by_source <- 
  reg_run(c("story_source", 
            "pres_approval_cat", "story_order"), 
          interaction_terms = c("story_source", 
                                "pres_approval_cat"),
          story_dummies = T,
          contrast_var = "story_source",
          quantities_calculate = "marginal means")

# estimates (marginal means) for the paper/Figure B5
# this is mean belief in stories, conditional on presidential approval
# and assigned news source
main_treatment_by_source_estimates <- main_treatment_by_source$effect |>
  # drop RBC and no-source stories: not considered in this analysis
  filter(!story_source %in% c("RBC", "No source")) |>
  mutate(story_source = factor(story_source, 
                               levels = c("Channel One", "Russia-24", 
                                          "RIA", "RT", "KP",
                                          "Meduza", "Rain", 
                                          "Echo of Moscow")),
         pres_approval_cat = case_when(
           pres_approval_cat == "Certainly disapprove" ~ "Strong\ncritic",
           pres_approval_cat == "Somewhat disapprove" ~ "Moderate\ncritic",
           pres_approval_cat == "Somewhat approve" ~ "Moderate\nsupporter",
           pres_approval_cat == "Certainly approve" ~ "Strong\nsupporter"
         ),
         pres_approval_cat = factor(pres_approval_cat, 
                                    levels = c("Strong\ncritic",
                                               "Moderate\ncritic",
                                               "Moderate\nsupporter",
                                               "Strong\nsupporter")),
         across(c(estimate, conf.low, conf.high), ~100*.x))

write_csv(main_treatment_by_source_estimates,
          "Estimates/Figure_B5_main_treatment_by_source_estimates.csv")

# Figure B5
ggplot(main_treatment_by_source_estimates, 
       aes(pres_approval_cat, estimate,
           color = story_source, shape = story_source)) +
  geom_point(size = 2,
             position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(ymin = conf.low,
                    ymax = conf.high),
                width = 0,
                position = position_dodge(width = 0.4)) +
  scale_color_manual(name = "Source", values = c("#a50026",
                                                 "#d73027",
                                                 "#f46d43",
                                                 "#fdae61",
                                                 "#fee090",
                                                 "#313695",
                                                 "#4575b4",
                                                 "#74add1")) +
  scale_shape_manual(name = "Source", values = c(17, 17, 17, 17, 17, 19, 19, 19)) +
  labs(x = "", y = "Probability of believing news stories, percent",
       title = "Perceived truthfulness of messages given news source") +
  ylim(c(35, 65)) +
  coord_flip() + 
  theme_classic() +
  theme(legend.position = "bottom")
ggsave("Figures/Figure_B5.png")

# Table B5
texreg(main_treatment_by_source$model,
       digits = 3,
       omit.coef = "survey_day|story_label",
       label = "tabonlinesourceeffect",
       custom.coef.names = c("Intercept", "Source: Meduza", "Source: Rain",
                             "Source: Echo of Moscow", "Source: RBC", 
                             "Source: Channel One", "Source: Russia-24", 
                             "Source: RT", "Source: RIA", "Source: KP",
                             "Somewhat disapprove", "Somewhat approve",
                             "Certainly approve",
                             "Story order", 
                             "Source: Meduza*Somewhat disapprove", 
                             "Source: Rain*Somewhat disapprove",
                             "Source: Echo of Moscow*Somewhat disapprove", 
                             "Source: RBC*Somewhat disapprove", 
                             "Source: Channel One*Somewhat disapprove", 
                             "Source: Russia-24*Somewhat disapprove", 
                             "Source: RT*Somewhat disapprove", 
                             "Source: RIA*Somewhat disapprove", 
                             "Source: KP*Somewhat disapprove", 
                             "Source: Meduza*Somewhat approve", 
                             "Source: Rain*Somewhat approve",
                             "Source: Echo of Moscow*Somewhat approve", 
                             "Source: RBC*Somewhat approve", 
                             "Source: Channel One*Somewhat approve", 
                             "Source: Russia-24*Somewhat approve", 
                             "Source: RT*Somewhat approve", 
                             "Source: RIA*Somewhat approve", 
                             "Source: KP*Somewhat approve", 
                             "Source: Meduza*Certainly approve", 
                             "Source: Rain*Certainly approve",
                             "Source: Echo of Moscow*Certainly approve", 
                             "Source: RBC*Certainly approve", 
                             "Source: Channel One*Certainly approve", 
                             "Source: Russia-24*Certainly approve", 
                             "Source: RT*Certainly approve", 
                             "Source: RIA*Certainly approve", 
                             "Source: KP*Certainly approve"),
       caption = paste("Treatment effect in the main study",
                       "(individual news sources)"),
       custom.note = paste("\\item %stars. Estimates from linear regressions",
                           "with news story evaluations as dependent variables.",
                           "Data from the social media sample.",
                           "The reference category in presidential approval",
                           "is 'Certainly disapprove.'",
                           "The reference category in source treatments",
                           "is 'No source.'",
                           "Story and day fixed effects included.",
                           "Heteroskedasticity-robust standard errors",
                           "in parentheses (clustered on respondent)."),
       float.pos = "H",
       fontsize = "scriptsize",
       single.row = T,
       caption.above = TRUE,
       threeparttable = TRUE,
       use.packages = FALSE,
       file = "Tables/Table_B5_main_treatment_by_source.tex",
       include.ci = F)

# APPENDIX Figure B6 and Table B8 ("Model 2")  ----------
# Experimental Results for Pre-Selected and “Current” News Stories

# a version of the model in Figure 3 where the effect is conditional
# on whether the story is pre-selected or current
main_treatment_effect_by_current <- 
  reg_run(c("state_controlled", 
            "pres_approval_cat", "story_order", "story_current"), 
          interaction_terms = c("state_controlled", 
                                "pres_approval_cat", 
                                "story_current"),
          story_dummies = F)

# estimates for the paper/Figure B6
main_treatment_effect_by_current_est <- 
  main_treatment_effect_by_current$effect |>
  mutate(pres_approval_cat = case_when(
    pres_approval_cat == "Certainly disapprove" ~ "Strong\ncritic",
    pres_approval_cat == "Somewhat disapprove" ~ "Moderate\ncritic",
    pres_approval_cat == "Somewhat approve" ~ "Moderate\nsupporter",
    pres_approval_cat == "Certainly approve" ~ "Strong\nsupporter"
  ),
  pres_approval_cat = factor(pres_approval_cat, 
                             levels = c("Strong\ncritic",
                                        "Moderate\ncritic",
                                        "Moderate\nsupporter",
                                        "Strong\nsupporter")),
  across(c(estimate, conf.low, conf.high), ~100*.x))

# save for the paper/Figure B6
write_csv(main_treatment_effect_by_current_est,
          "Estimates/Figure_B6_main_treatment_effect_current_estimates.csv")

# Figure B6
ggplot(main_treatment_effect_by_current_est, 
       aes(pres_approval_cat, estimate,
           shape = story_current, linetype = story_current)) +
  geom_point(size = 2,
             position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(ymin = conf.low,
                    ymax = conf.high),
                width = 0,
                position = position_dodge(width = 0.4),
                size = 0.8) +
  geom_hline(yintercept = 0, linetype = "dashed", col = "gray") +
  labs(x = "", y = "Effect of state vs independent media",
       title = "The effect of state media by story type",
       linetype = "Stories", shape = "Stories") +
  ylim(c(-20, 20)) +
  theme_classic() +
  theme(legend.position = "bottom") +
  coord_flip()

ggsave("Figures/Figure_B6.png")

# regression Table B8
# (the model by story direction/Model 1 was estimated earlier, for Figure 3)
texreg(list("Model 1" = main_treatment_effect_by_direction$model,
            "Model 2" = main_treatment_effect_by_current$model),
       digits = 3,
       omit.coef = "survey_day|story_label",
       custom.coef.names = c("Intercept", "Source: Critical", 
                             "Source: State-controlled",
                             "Source: RBC", "Somewhat disapprove",
                             "Somewhat approve", "Certainly approve",
                             "Story order", "Neutral story",
                             "Pro-Russia story",
                             "Source: Critical*Somewhat disapprove",
                             "Source: State-controlled*Somewhat disapprove",
                             "Source: RBC*Somewhat disapprove",
                             "Source: Critical*Somewhat approve",
                             "Source: State-controlled*Somewhat approve",
                             "Source: RBC*Somewhat approve",
                             "Source: Critical*Certainly approve",
                             "Source: State-controlled*Certainly approve",
                             "Source: RBC*Certainly approve",
                             "Source: Critical*Neutral story",
                             "Source: State-controlled*Neutral story",
                             "Source: RBC*Neutral story",
                             "Source: Critical*Pro-Russia story",
                             "Source: State-controlled*Pro-Russia story",
                             "Source: RBC*Pro-Russia story",
                             "Somewhat disapprove*Neutral story",
                             "Somewhat approve*Neutral story",
                             "Certainly approve*Neutral story",
                             "Somewhat disapprove*Pro-Russia story",
                             "Somewhat approve*Pro-Russia story",
                             "Certainly approve*Pro-Russia story",
                             "Source: Critical*Somewhat disapprove*Neutral story",
                             "Source: State-controlled*Somewhat disapprove*Neutral story",
                             "Source: RBC*Somewhat disapprove*Neutral story",
                             "Source: Critical*Somewhat approve*Neutral story",
                             "Source: State-controlled*Somewhat approve*Neutral story",
                             "Source: RBC*Somewhat approve*Neutral story",
                             "Source: Critical*Certainly approve*Neutral story",
                             "Source: State-controlled*Certainly approve*Neutral story",
                             "Source: RBC*Certainly approve*Neutral story",
                             "Source: Critical*Somewhat disapprove*Pro-Russia story",
                             "Source: State-controlled*Somewhat disapprove*Pro-Russia story",
                             "Source: RBC*Somewhat disapprove*Pro-Russia story",
                             "Source: Critical*Somewhat approve*Pro-Russia story",
                             "Source: State-controlled*Somewhat approve*Pro-Russia story",
                             "Source: RBC*Somewhat approve*Pro-Russia story",
                             "Source: Critical*Certainly approve*Pro-Russia story",
                             "Source: State-controlled*Certainly approve*Pro-Russia story",
                             "Source: RBC*Certainly approve*Pro-Russia story",
                             "Pre-selected story",
                             "Source: Critical*Pre-selected story",
                             "Source: State-controlled*Pre-selected story",
                             "Source: RBC*Pre-selected story",
                             "Somewhat disapprove*Pre-selected story",
                             "Somewhat approve*Pre-selected story",
                             "Certainly approve*Pre-selected story",
                             "Source: Critical*Somewhat disapprove*Pre-selected story",
                             "Source: State-controlled*Somewhat disapprove*Pre-selected story",
                             "Source: RBC*Somewhat disapprove*Pre-selected story",
                             "Source: Critical*Somewhat approve*Pre-selected story",
                             "Source: State-controlled*Somewhat approve*Pre-selected story",
                             "Source: RBC*Somewhat approve*Pre-selected story",
                             "Source: Critical*Certainly approve*Pre-selected story",
                             "Source: State-controlled*Certainly approve*Pre-selected story",
                             "Source: RBC*Certainly approve*Pre-selected story"),
       label = "tabonlinestatemediaeffecttopic",
       caption = paste("Treatment effect in the main study",
                       "given news story content"),
       custom.note = paste("\\item %stars. Estimates from linear regressions",
                           "with news story evaluations as dependent variables.",
                           "Data from the social media sample.",
                           "The reference category in presidential approval",
                           "is 'Certainly disapprove.'",
                           "The reference category in story content in Model 1",
                           "is 'Critical story.'",
                           "The reference category in story content in Model 2",
                           "is 'Current story.'",
                           "Day fixed effects included.",
                           "Heteroskedasticity-robust standard errors",
                           "in parentheses (clustered on respondent)."),
       caption.above = TRUE,
       float.pos = "H",
       fontsize = "tiny",
       threeparttable = TRUE,
       use.packages = FALSE,
       single.row = T,
       file = "Tables/Table_B8_main_treatment_effect_topic.tex",
       include.ci = F)

# APPENDIX Figure B7 and Table B6  ---------------- 
# Alternative Categorizations of State-Controlled Media Outlets 
# model from Figure 3 with an alternative definition of source treatment
# where RBC is considered a state-controlled outlet
main_treatment_effect_rbc <- 
  reg_run(c("state_controlled_RBC",
            "pres_approval_cat", "story_order"), 
          interaction_terms = c("state_controlled_RBC", 
                                "pres_approval_cat"),
          story_dummies = T,
          contrast_var = "state_controlled_RBC")

# model from Figure 3 with an alternative definition of source treatment
# where only state-run (state-owned) outlets are considered state-controlled
main_treatment_effect_state_run <- 
  reg_run(c("state_run",
            "pres_approval_cat", "story_order"), 
          interaction_terms = c("state_run", 
                                "pres_approval_cat"),
          story_dummies = T,
          contrast_var = "state_run",
          contrast_filter = "Critical - (State-run)")

# estimates for Figure B7
main_treatment_effect_alt_estimates <- bind_rows(
  main_treatment_effect_rbc$effect |>
    mutate(state_definition = "RBC as a state media outlet"),
  main_treatment_effect_state_run$effect |>
    mutate(state_definition = "Only state-owned outlets")) |>
  mutate(pres_approval_cat = case_when(
    pres_approval_cat == "Certainly disapprove" ~ "Strong\ncritic",
    pres_approval_cat == "Somewhat disapprove" ~ "Moderate\ncritic",
    pres_approval_cat == "Somewhat approve" ~ "Moderate\nsupporter",
    pres_approval_cat == "Certainly approve" ~ "Strong\nsupporter"
  ),
  pres_approval_cat = factor(pres_approval_cat, 
                             levels = c("Strong\ncritic",
                                        "Moderate\ncritic",
                                        "Moderate\nsupporter",
                                        "Strong\nsupporter")),
  across(c(estimate, conf.low, conf.high), ~100*.x))

# save estimates for the paper
write_csv(main_treatment_effect_alt_estimates,
          "Estimates/Figure_B7_main_treatment_effect_alt_def_estimates.csv")

# Figure B7
ggplot(main_treatment_effect_alt_estimates, 
       aes(pres_approval_cat, estimate, 
           linetype = state_definition, shape = state_definition)) +
  geom_point(position = position_dodge(width = 0.4), size = 2) +
  geom_errorbar(aes(ymin = conf.low,
                    ymax = conf.high),
                width = 0,
                position = position_dodge(width = 0.4),
                size = 0.8) +
  geom_hline(yintercept = 0, linetype = "dashed", col = "gray") +
  labs(x = "", y = "Effect of state vs independent media",
       title = "The effect of state media",
       subtitle = "For alternative definitions of state-controlled media",
       linetype = "", shape = "") +
  ylim(c(-20, 20)) +
  coord_flip() +
  theme_classic() +
  theme(legend.position = "bottom")

# Table B6
texreg(list("Model 1" = main_treatment_effect_rbc$model,
            "Model 2" = main_treatment_effect_state_run$model),
       digits = 3,
       omit.coef = "survey_day|story_label",
       custom.coef.names = c("Intercept", "Source: Critical", 
                             "Source: State-controlled",
                             "Somewhat disapprove", "Somewhat approve",
                             "Certainly approve", "Story order", 
                             "Source: Critical*Somewhat disapprove",
                             "Source: State-controlled*Somewhat disapprove",
                             "Source: Critical*Somewhat approve",
                             "Source: State-controlled*Somewhat approve",
                             "Source: Critical*Certainly approve",
                             "Source: State-controlled*Certainly approve",
                             "Source: Critical",
                             "Source: State-owned",
                             "Source: Other",
                             "Source: Critical*Somewhat disapprove",
                             "Source: State-owned*Somewhat disapprove",
                             "Source: Other*Somewhat disapprove",
                             "Source: Critical*Somewhat approve",
                             "Source: State-owned*Somewhat approve",
                             "Source: Other*Somewhat approve",
                             "Source: Critical*Certainly approve",
                             "Source: State-owned*Certainly approve",
                             "Source: Other*Certainly approve"),
       label = "tabonlinestatemediaeffectalt",
       caption = paste("Treatment effect in the main study",
                       "(alternative definitions of state-controlled media)"),
       custom.note = paste("\\item %stars. Estimates from linear regression models.",
                           "with news story evaluations as dependent variables.",
                           "In Model 1, RBC is treated as a state-controlled outlet.",
                           "In Model 2, state-controlled outlets are divided",
                           "into 'State-owned' and 'Other.'",
                           "In Model 1, RBC is treated as a state-controlled outlet.",
                           "The reference category in presidential approval",
                           "is 'Certainly disapprove.'",
                           "The reference category in source treatments",
                           "is 'No source.'",
                           "Data from the social media sample.",
                           "Story and day fixed effects included.",
                           "Heteroskedasticity-robust standard errors",
                           "in parentheses (clustered on respondent)."),
       #       single.row = T,
       fontsize = "scriptsize",
       float.pos = "H",
       caption.above = TRUE,
       threeparttable = TRUE,
       use.packages = FALSE,
       file = "Tables/Table_B6_main_treatment_effect_alt.tex",
       include.ci = FALSE)


# APPENDIX Table C1 and Figure C1 ----------
# regressions based on the Levada survey

# treatment effect conditional on Putin approval
levada_treatment_effect <- lm_robust(
  story_response_scaled ~ story_source_tv1 +
    female + age  +
    education + story_label + 
    pres_approval_cat*story_source_tv1,
  data = levada_stories,
  weights = weight,
  clusters = respondent_id,
  se_type = "CR0")

# treatment effect conditional on vote choice in 2018
levada_treatment_effect_vote <- lm_robust(
  story_response_scaled ~ story_source_tv1 +
    female + age +
    education + story_label + 
    vote_outcome_2018*story_source_tv1,
  data = levada_stories,
  weights = weight,
  clusters = respondent_id,
  se_type = "CR0")

# combine estimates (contrasts) for paper/Figure C1
levada_treatment_effect_est <- emmeans(
  levada_treatment_effect,
  ~pres_approval_cat*story_source_tv1)
# contrasts for the figure
levada_treatment_effect_est <- summary(pairs(
  levada_treatment_effect_est, simple = "story_source_tv1")) |>
  # reverse the sign so that the contrast is "state media" - "independent media"
  mutate(estimate = -estimate,
         conf.low = estimate - 1.96*SE,
         conf.high = estimate + 1.96*SE,
         pres_approval_cat = case_when(
           pres_approval_cat == "Certainly disapprove" ~ "Strong\ncritic",
           pres_approval_cat == "Somewhat disapprove" ~ "Moderate\ncritic",
           pres_approval_cat == "Somewhat approve" ~ "Moderate\nsupporter",
           pres_approval_cat == "Certainly approve" ~ "Strong\nsupporter"
         ),
         pres_approval_cat = factor(pres_approval_cat, 
                                    levels = c("Strong\ncritic",
                                               "Moderate\ncritic",
                                               "Moderate\nsupporter",
                                               "Strong\nsupporter")),
         across(c(estimate, conf.low, conf.high), ~100*.x))

# contrasts for the figure
levada_treatment_effect_vote_est <- emmeans(
  levada_treatment_effect_vote,
  ~vote_outcome_2018*story_source_tv1)
levada_treatment_effect_vote_est <- summary(pairs(
  levada_treatment_effect_vote_est, simple = "story_source_tv1")) |>
  # reverse the sign so that the contrast is "state media" - "independent media"
  mutate(estimate = -estimate,
         conf.low = estimate - 1.96*SE,
         conf.high = estimate + 1.96*SE,
         highlight = case_when(
           vote_outcome_2018 %in% c("Putin", "Liberal") ~ "black",
           TRUE ~ "gray"
         ),
         vote_outcome_2018 = case_when(
           vote_outcome_2018 == "Putin" ~ "Putin supporter",
           vote_outcome_2018 == "Liberal" ~ "Liberal critic",
           vote_outcome_2018 == "Communist/nationalist" ~ 
             "Illiberal critic\n(communist,\nnationalist)",
           TRUE ~ vote_outcome_2018
         ),
         across(c(estimate, conf.low, conf.high), ~100*.x))

# save estimates for the paper
write_csv(levada_treatment_effect_est, 
          "Estimates/Figure_C1_levada_treatment_effect_estimates.csv")
write_csv(levada_treatment_effect_vote_est, 
          "Estimates/Figure_C1_levada_treatment_effect_vote_estimates.csv")

# Figure C1
c1left <- ggplot(levada_treatment_effect_vote_est, 
                 aes(vote_outcome_2018, estimate,
                     color = highlight)) +
  geom_point(size = 2,
             position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(ymin = conf.low,
                    ymax = conf.high),
                width = 0,
                position = position_dodge(width = 0.4),
                size = 0.6) +
  scale_color_manual(values = c("black", "gray")) +
  geom_hline(yintercept = 0, linetype = "dashed", col = "gray") +
  labs(x = "", y = "Effect of state media",
       title = "Given vote choice in 2018") +
  ylim(c(-27, 27)) +
  theme_classic() +
  theme(legend.position = "none") +
  coord_flip()

c1right <- ggplot(levada_treatment_effect_est, 
                  aes(pres_approval_cat, estimate)) +
  geom_point(size = 2,
             position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(ymin = conf.low,
                    ymax = conf.high),
                width = 0,
                position = position_dodge(width = 0.4),
                size = 0.6) +
  geom_hline(yintercept = 0, linetype = "dashed", col = "gray") +
  labs(x = "", y = "Effect of state media",
       title = "Given Putin approval") +
  ylim(c(-27, 27)) +
  theme_classic() +
  coord_flip()

c1left + c1right +
  plot_annotation(title = "The effect of state media in the national sample")
ggsave("Figures/Figure_C1.png")

# Table C1
texreg(list("Model 1" = levada_treatment_effect,
            "Model 2" = levada_treatment_effect_vote),
       digits = 3,
       float.pos = "H",
       custom.coef.names = c("Intercept", "Channel One", "Female",
                             "Age", "Education", 
                             "Somewhat disapprove", "Somewhat approve",
                             "Certainly approve", 
                             "Channel One*Somewhat disapprove",
                             "Channel One*Somewhat approve",
                             "Channel One*Certainly approve",
                             "Voted liberal",
                             "Voted for Putin",
                             "Spoiled ballot/no vote",
                             "Channel One*Liberal",
                             "Channel One*Putin",
                             "Channel One*No vote"),
       omit.coef = "story_label",
       label = "tablevadastatemediaeffect",
       caption = paste("Treatment effect in the nationally representative survey"),
       custom.note = paste("\\item %stars. Estimates from linear regressions",
                           "with news story evaluations as dependent variables.",
                           "In Model 1, regime support is measured",
                           "via presidential approval. In Model 2,",
                           "regime support is measured via vote outcome",
                           "in the 2018 presidential election.",
                           "The reference category in presidential approval",
                           "is 'Certainly disapprove.'",
                           "The reference category in 2018 vote",
                           "is 'Communist/nationalist.'",
                           "Data from the Levada sample.",
                           "Story fixed effects included.",
                           "Heteroskedasticity-robust standard errors",
                           "in parentheses (clustered on respondent)."),
       fontsize = "scriptsize",
       caption.above = TRUE,
       threeparttable = TRUE,
       use.packages = FALSE,
       file = "Tables/Table_C1_levada_treatment_effect.tex",
       include.ci = F)

# 3.4 Additional Evidence ---------------
# results from the OMI survey, in some cases in comparison to the main study
# and the Levada survey

# Figure 4 and Table D2 ------
# Putin supporters list state TV among trusted sources

# regressions of trust in state-controlled media/state TV/independent media
# on presidential approval and covariates
state_source_trust_omi <- media_attitudes("source_trusted_state_controlled", 
                                          "State-controlled",
                                          study = "omi")
state_tv_source_trust_omi <- media_attitudes("source_trusted_state_tv", 
                                             "State TV",
                                             study = "omi")
ind_source_trust_omi <- media_attitudes("source_trusted_independent", 
                                        "Independent",
                                        study = "omi")

# save estimates (conditional means) for the paper/Figure 4
source_trust_estimates <- bind_rows(state_tv_source_trust_omi$means,
                                    ind_source_trust_omi$means) |>
  mutate(x = case_when(
    x == "Certainly disapprove" ~ "Strong\ncritic",
    x == "Somewhat disapprove" ~ "Moderate\ncritic",
    x == "Somewhat approve" ~ "Moderate\nsupporter",
    x == "Certainly approve" ~ "Strong\nsupporter"
  ),
  x = factor(x, levels = c("Strong\ncritic",
                           "Moderate\ncritic",
                           "Moderate\nsupporter",
                           "Strong\nsupporter")),
  media = factor(media, levels = c("State TV", "Independent")),
  across(c(estimate, conf.low, conf.high), ~100*.x))

# save estimates for the paper/Figure 4
write_csv(source_trust_estimates,
          "Estimates/Figure_4_OMI_media_trust.csv")

# Figure 4
ggplot(source_trust_estimates, 
       aes(x, estimate, linetype = media, shape = media)) +
  geom_point(size = 2,
             position = position_dodge(width = 0.4)) +
  geom_text(aes(label = media, vjust = if_else(media == "State TV", -.5, -2.5)), 
            size = 3,
            show.legend = F) +
  geom_errorbar(aes(ymin = conf.low,
                    ymax = conf.high),
                width = 0,
                position = position_dodge(width = 0.4), size = 0.8) +
  labs(x = "", y = "Probability of trusting state TV or independent media, percent",
       title = "Putin supporters list state TV among trusted sources",
       linetype = "", shape = "") +
  ylim(c(-2, 90)) +
  theme_classic() +
  theme(legend.position = "bottom") +
  coord_flip()
ggsave("Figures/Figure_4.png")

# Table D2
# list of models
model_list_trust <- list("State media" = state_source_trust_omi$model,
                         "State TV" = state_tv_source_trust_omi$model,
                         "Independent" = ind_source_trust_omi$model)
# list of vcov matrices
vcov_list_trust <- list(state_source_trust_omi$vcov_hc,
                        state_tv_source_trust_omi$vcov_hc,
                        ind_source_trust_omi$vcov_hc)
# calculate standard errors based on above (for the table)
se_list_trust <- map2(model_list_trust, vcov_list_trust,
                      ~coeftest(.x, vcov. = .y)[, 2])
# calculate p-values based on above (for the table)
pval_list_trust <- map2(model_list_trust, vcov_list_trust,
                        ~coeftest(.x, vcov. = .y)[, 4])

texreg(model_list_trust,
       digits = 3,
       custom.coef.names = c("Intercept", "Somewhat disapprove",
                             "Somewhat approve", "Certainly approve",
                             "Age", "Female", "Education"),
       label = "tabomimediatrust",
       caption = paste("Trust in state and independent media"),
       override.se = se_list_trust,
       override.pvalues = pval_list_trust,
       custom.note = paste("\\item %stars. Estimates from a linear regression",
                           "with trust in state and independent media",
                           "as dependent variables.",
                           "Data from the OMI survey (media perceptions survey).",
                           "Heteroskedasticity-robust standard errors",
                           "in parentheses."),
       float.pos = "H",
       fontsize = "scriptsize",
       caption.above = TRUE,
       threeparttable = TRUE,
       use.packages = FALSE,
       file = "Tables/Table_D2_OMI_media_trust.tex")

# Figure 5 and Tables D3-D6 -------
# Supporters evaluate state media as accurate despite bias

# recode the data set for multinomial model analysis
omi_rankings <- omi_resp |>
  mutate(pres_approval_dummy = case_when(
    pres_approval_dummy == 0 ~ "Putin critic",
    pres_approval_dummy == 1 ~ "Putin supporter"
  ),
  pres_approval_dummy = factor(pres_approval_dummy,
                               levels = c("Putin supporter",
                                          "Putin critic"))) |>
  # keep only variables with media evaluations and necessary covariates
  select(starts_with("ranking_"), "pres_approval_cat",
         "pres_approval_dummy", "female", "age_ordinal", "education_higher",
         "source_trusted_independent") |>
  # drop the evaluations of Echo of Moscow and RBC (not considered in this analysis)
  select(-ends_with("Echo"), -ends_with("RBC")) |>
  mutate_at(vars(starts_with("ranking_")),
            list(~as.factor(.))) |>
  mutate_at(vars(starts_with("ranking_")),
            list(~relevel(., ref = "Don't know/hard to say")))

# list of dependent variables for regression models
# (evaluations of accuracy, bias, etc., for the selected state media outlets)
ranking_list <- omi_rankings |> 
  select(starts_with("ranking")) |> 
  names()

# run multinomial models on each dependent variable from the list
# in each model, the dependent variable is the categorical variable
# evaluating accuracy, bias, or another characteristic of a particular outlet 
# (Channel One, RT, etc.)
rankings_regs <- map(
  ranking_list, 
  ~run_multinom(.x, reg_formula = " ~ pres_approval_dummy + female + 
                      age_ordinal + education_higher"))

# extract the list of models
rankings_mods <- map(rankings_regs, "model")

# extract conditional means/probabilities from objects created by the function
rankings_probs <- map_df(rankings_regs, "probs") |>
  # drop "hard to say" responses
  filter(response != "Don.t.know.hard.to.say") |>
  # create a variable indicating whether the evaluation is positive or negative
  # (a technical indicator used to select estimates for plots)
  mutate(eval_direction = case_when(
    response == "Mostly.yes" ~ "positive",
    response %in% c("Omit.important", "Often.gives.false.info", 
                    "Pro.government", "Anti.government", 
                    "Not.independent") ~ "negative",
  )) |>
  filter(eval_direction == "positive") |>
  mutate(across(c(prob, conf.low, conf.high), ~100*.x),
         prob = round(prob, 1),
         source = factor(source, 
                         levels = c("RT", "RIA", "Russia-24", 
                                    "Channel One")),
         dimension = case_when(
           dimension == "accuracy" ~ "Accurate",
           dimension == "independence" ~ "Politically\nindependent",
           dimension == "bias" ~ "Unbiased,\nneutral",
           dimension == "completeness" ~ "Not censored"
         ))

# save estimates for the paper/Figure 5
write_csv(rankings_probs, "Estimates/Figure_5_omi_perceptions_probs_all.csv")

# Figure 5
ggplot(rankings_probs, 
       aes(source, prob, shape = pres_approval_dummy,
           linetype = pres_approval_dummy, color = pres_approval_dummy)) +
  geom_point(size = 1.5) +
  geom_errorbar(aes(ymin = conf.low,
                    ymax = conf.high),
                width = 0, size = 0.7) +
  geom_text(aes(label = prob), vjust = -1, size = 3, show.legend = F) +
  scale_color_manual(values = c("black", "darkgray")) +
  labs(x = "", 
       y = "Probability of saying that media outlets are\naccurate, uncensored, independent, or politically unbiased",
       title = "Supporters evaluate state media as accurate despite bias",
       linetype = "", shape = "", color = "") +
  theme_classic() +
  coord_flip() +
  theme(legend.position = "bottom") +
  facet_wrap(~dimension, ncol = 4)
ggsave("Figures/Figure_5.png")

# regression tables
# add news source labels to regression models
source_names <- c("RT", "Channel 1", "Russia-24", "RIA")
completeness_mods <- rankings_mods[1:4]
names(completeness_mods) <- source_names
accuracy_mods <- rankings_mods[5:8]
names(accuracy_mods) <- source_names
bias_mods <- rankings_mods[9:12]
names(bias_mods) <- source_names
independence_mods <- rankings_mods[13:16]
names(independence_mods) <- source_names

# Table D3
texreg(completeness_mods,
       digits = 3,
       custom.coef.names = c("Y: Intercept", 
                             "Y: Critic",
                             "Y: Female", 
                             "Y: Age", 
                             "Y: Higher education",
                             "N: Intercept", 
                             "N: Critic",
                             "N: Female", 
                             "N: Age", 
                             "N: Higher education"),
       label = "tabomirankingcompleteness",
       caption = paste("State media evaluations:",
                       "Completeness"),
       custom.note = paste("\\item %stars. Estimates from multinomial regressions",
                           "with media evaluations as dependent variables.",
                           "Outcomes: Mostly complete (Y in the table),", 
                           "Omits important information (N in the table),",
                           "and Hard to say (reference category).",
                           "Data from the OMI survey (media perceptions survey).",
                           "Standard errors in parentheses."),
       float.pos = "H",
       fontsize = "tiny",
       single.row = T,
       caption.above = TRUE,
       threeparttable = TRUE,
       use.packages = FALSE,
       file = "Tables/Table_D3_OMI_ranking_completeness.tex")

# Table D4
texreg(accuracy_mods,
       digits = 3,
       custom.coef.names = c("Y: Intercept", 
                             "Y: Critic",
                             "Y: Female", 
                             "Y: Age", 
                             "Y: Higher education",
                             "N: Intercept", 
                             "N: Critic",
                             "N: Female", 
                             "N: Age", 
                             "N: Higher education"),
       label = "tabomirankingaccuracy",
       caption = paste("State media evaluations:",
                       "Accuracy"),
       custom.note = paste("\\item %stars. Estimates from multinomial regressions",
                           "with media evaluations as dependent variables.",
                           "Outcomes: Mostly accurate (Y in the table),", 
                           "Often gives false information (N in the table),",
                           "and Hard to say (reference category).",
                           "Data from the OMI survey (media perceptions survey).",
                           "Standard errors in parentheses."),
       float.pos = "H",
       fontsize = "tiny",
       single.row = T,
       caption.above = TRUE,
       threeparttable = TRUE,
       use.packages = FALSE,
       file = "Tables/Table_D4_OMI_ranking_accuracy.tex")

# Table D5
texreg(independence_mods,
       digits = 3,
       custom.coef.names = c("Y: Intercept", 
                             "Y: Critic",
                             "Y: Female", 
                             "Y: Age", 
                             "Y: Higher education",
                             "N: Intercept", 
                             "N: Critic",
                             "N: Female", 
                             "N: Age", 
                             "N: Higher education"),
       label = "tabomirankingindependence",
       caption = paste("State media evaluations:",
                       "Independence"),
       custom.note = paste("\\item %stars. Estimates from multinomial regressions",
                           "with media evaluations as dependent variables.",
                           "Outcomes: Mostly independent from authorities",
                           "(Y in the table),", 
                           "Not independent (N in the table),",
                           "and Hard to say (reference category).",
                           "Data from the OMI survey (media perceptions survey).",
                           "Standard errors in parentheses."),
       float.pos = "H",
       fontsize = "tiny",
       single.row = T,
       caption.above = TRUE,
       threeparttable = TRUE,
       use.packages = FALSE,
       file = "Tables/Table_D5_OMI_ranking_independence.tex")

# Table D6
texreg(bias_mods,
       digits = 3,
       custom.coef.names = c("Anti: Intercept", 
                             "Anti: Critic",
                             "Anti: Female", 
                             "Anti: Age", 
                             "Anti: Higher education",
                             "Y: Intercept", 
                             "Y: Critic",
                             "Y: Female", 
                             "Y: Age", 
                             "Y: Higher education",
                             "Pro: Intercept", 
                             "Pro: Critic",
                             "Pro: Female", 
                             "Pro: Age", 
                             "Pro: Higher education"),
       label = "tabomirankingbias",
       caption = paste("State media evaluations:",
                       "Political bias"),
       custom.note = paste("\\item %stars. Estimates from multinomial regressions",
                           "with media evaluations as dependent variables.",
                           "Outcomes: Mostly neutral (Y in the table),", 
                           "Anti-government (Anti in the table),",
                           "Pro-government (Pro in the table),",
                           "and Hard to say (reference category).",
                           "Data from the OMI survey (media perceptions survey).",
                           "Standard errors in parentheses."),
       float.pos = "H",
       fontsize = "tiny",
       single.row = T,
       caption.above = TRUE,
       threeparttable = TRUE,
       use.packages = FALSE,
       file = "Tables/Table_D6_OMI_ranking_bias.tex")


# APPENDIX Figure D1 and Table D1 -------
# State and independent media usage

# regressions of dummies indicating usage of state-controlled outlets,
# state TV, or independent media, on presidential approval and other covariates

# main study
state_source_usage <- media_attitudes("source_state_controlled_dummy", 
                                      "State-controlled",
                                      study = "main")
state_tv_source_usage <- media_attitudes("source_state_tv", 
                                         "State TV",
                                         study = "main")
ind_source_usage <- media_attitudes("source_independent_dummy", 
                                    "Foreign/critical media",
                                    study = "main")
# OMI
state_source_usage_omi <- media_attitudes("source_used_state_controlled", 
                                          "State-controlled",
                                          study = "omi")
state_tv_source_usage_omi <- media_attitudes("source_used_state_tv", 
                                             "State TV",
                                             study = "omi")
ind_source_usage_omi <- media_attitudes("source_used_independent", 
                                        "Foreign/critical media",
                                        study = "omi")
# Levada
state_tv_source_usage_levada <- media_attitudes("info_federal_state_tv_dummy", 
                                                "State TV",
                                                study = "levada")
ind_source_usage_levada <- media_attitudes("info_indep_media_dummy", 
                                           "Foreign/critical media",
                                           study = "levada")

# estimates for the paper/Figure D1
state_indep_usage <- bind_rows(state_tv_source_usage_levada$means |>
                                 mutate(survey = "National"),
                               ind_source_usage_levada$means |>
                                 mutate(survey = "National"),
                               state_tv_source_usage_omi$means |>
                                 mutate(survey = "Online panel"),
                               ind_source_usage_omi$means |>
                                 mutate(survey = "Online panel"),
                               state_tv_source_usage$means |>
                                 mutate(survey = "Main study"),
                               ind_source_usage$means |>
                                 mutate(survey = "Main study")) |>
  select(-group) |>
  mutate(x = case_when(
    x == "Certainly disapprove" ~ "Strong\ncritic",
    x == "Somewhat disapprove" ~ "Moderate\ncritic",
    x == "Somewhat approve" ~ "Moderate\nsupporter",
    x == "Certainly approve" ~ "Strong\nsupporter"
  ),
  x = factor(x, levels = c("Strong\ncritic",
                           "Moderate\ncritic",
                           "Moderate\nsupporter",
                           "Strong\nsupporter")))

# save estimates for the paper
write_csv(state_indep_usage, "Estimates/Figure_D1_media_usage_3_surveys.csv")

# Figure D1
ggplot(state_indep_usage, 
       aes(x, estimate, linetype = survey, shape = survey)) +
  geom_point(size = 2,
             position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(ymin = conf.low,
                    ymax = conf.high),
                width = 0,
                position = position_dodge(width = 0.4),
                size = 0.6) +
  labs(x = "", y = "Probability of using state TV or independent media",
       title = "Usage of state and independent media",
       linetype = "Sample", shape = "Sample") +
  ylim(c(0, 1)) +
  coord_flip() +
  theme_classic() +
  theme(legend.position = "bottom") +
  facet_wrap(~media, ncol = 2)
ggsave("Figures/Figure_D1.png")

# Table D1
# list of models and vcov matrices
model_list_usage <- list("State media" = state_source_usage$model,
                         "State TV" = state_tv_source_usage$model,
                         "Independent" = ind_source_usage$model,
                         "State media" = state_source_usage_omi$model,
                         "State TV" = state_tv_source_usage_omi$model,
                         "Independent" = ind_source_usage_omi$model,
                         "State TV" = state_tv_source_usage_levada$model,
                         "Independent" = ind_source_usage_levada$model)
vcov_list_usage <- list(state_source_usage$vcov_hc,
                        state_tv_source_usage$vcov_hc,
                        ind_source_usage$vcov_hc,
                        state_source_usage_omi$vcov_hc,
                        state_tv_source_usage_omi$vcov_hc,
                        ind_source_usage_omi$vcov_hc,
                        state_tv_source_usage_levada$vcov_hc,
                        ind_source_usage_levada$vcov_hc)
# calculate robust standard errors and p-values for the table
se_list_usage <- map2(model_list_usage, vcov_list_usage,
                      ~coeftest(.x, vcov. = .y)[, 2])
pval_list_usage <- map2(model_list_usage, vcov_list_usage,
                        ~coeftest(.x, vcov. = .y)[, 4])

texreg(model_list_usage,
       digits = 3,
       custom.coef.names = c("Intercept", "Somewhat disapprove",
                             "Somewhat approve", "Certainly approve",
                             "Age", "Female", "Education"),
       custom.header = list("Main study" = 1:3, "OMI Survey" = 4:6, 
                            "National survey" = 7:8),
       label = "tabmediausage",
       caption = paste("State and independent media usage"),
       override.se = se_list_usage,
       override.pvalues = pval_list_usage,
       custom.note = paste("\\item %stars. Estimates from a linear regression",
                           "with state and independent media usage",
                           "as dependent variables.",
                           "Data from the main study, media perceptions survey",
                           "(OMI), and the nationally representative survey",
                           "(Levada). In the regressions for the main study",
                           "and for the OMI survey, education is",
                           "dichotomized, and age is an ordinal measure.",
                           "Heteroskedasticity-robust standard errors",
                           "in parentheses."),
       float.pos = "H",
       fontsize = "scriptsize",
       caption.above = TRUE,
       threeparttable = TRUE,
       use.packages = FALSE,
       file = "Tables/Table_D1_media_usage_3_surveys.tex")

# APPENDIX Figure D2 ---------------
# Knowledge of Independent Media and Trust in/Usage of State Media

# regression of the dummy indicating the usage of state TV 
# conditional on presidential approval and knowledge of independent media
state_tv_usage_omi_aw <- lm(
  source_used_state_tv ~
    pres_approval_cat + female + age_ordinal + education_higher +
    pres_approval_cat*source_known_independent,
  data = omi_resp |>
    drop_na(pres_approval_cat, source_known_independent))

# regression of the dummy indicating trust in state TV 
# conditional on presidential approval and knowledge of independent media
state_tv_trust_omi_aw <- lm(
  source_trusted_state_tv ~
    pres_approval_cat + female + age_ordinal + education_higher +
    pres_approval_cat*source_known_independent,
  data = omi_resp |>
    drop_na(pres_approval_cat, source_known_independent))

# marginal means for the figure (trust conditional on approval and knowledge
# of independent media)
# usage of state TV
state_tv_usage_omi_aw_mm <- ggpredict(
  state_tv_usage_omi_aw,
  terms = c("pres_approval_cat", "source_known_independent"),
  vcov.fun = "vcovHC", vcov.type = "HC0") |>
  filter(x == "Certainly approve" | x == "Somewhat approve") 

# trust in state TV
state_tv_trust_omi_aw_mm <- ggpredict(
  state_tv_trust_omi_aw,
  terms = c("pres_approval_cat", "source_known_independent"),
  vcov.fun = "vcovHC", vcov.type = "HC0") |>
  filter(x == "Certainly approve" | x == "Somewhat approve")

# clean up estimates for the plot
state_tv_omi_aw_mm <- bind_rows(
  as.data.frame(state_tv_usage_omi_aw_mm) |> 
    mutate(outcome = "Usage of state TV"),
  as.data.frame(state_tv_trust_omi_aw_mm) |> 
    mutate(outcome = "Trust in state TV")
) |>
  mutate(knowledge_indep = case_when(
    group == 1 ~ "Yes",
    group == 0 ~ "No"
  ),
  x = case_when(
    x == "Somewhat approve" ~ "Moderate\nsupporter",
    x == "Certainly approve" ~ "Strong\nsupporter"
  ),
  x = factor(x, levels = c("Moderate\nsupporter",
                           "Strong\nsupporter")),
  knowledge_indep = factor(knowledge_indep, 
                           levels = c("Yes", 
                                      "No"))) |>
  rename(estimate = predicted)

# save estimates for the paper
write_csv(state_tv_omi_aw_mm,
          "Estimates/Figure_D2_omi_state_tv_by_awareness.csv")

# Figure D2
ggplot(state_tv_omi_aw_mm, 
       aes(x, estimate, linetype = knowledge_indep, 
           shape = knowledge_indep)) +
  geom_point(size = 2,
             position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(ymin = conf.low,
                    ymax = conf.high),
                width = 0,
                position = position_dodge(width = 0.4), size = 0.7) +
  labs(x = "", y = "Probability of trusting or using state TV",
       title = "Trust in/usage of state TV among supporters",
       subtitle = "Given knowledge of independent media",
       linetype = "Knows some independent media", 
       shape = "Knows some independent media") +
  ylim(c(0, 1)) +
  theme_classic() +
  theme(legend.position = "bottom") +
  coord_flip() +
  facet_wrap(~outcome, ncol = 2)
ggsave("Figures/Figure_D2.png")

# APPENDIX Figure D3 ----------------
# Knowledge of Independent Media and the Evaluations of State Media

# recode the data set for multinomial model analysis (same as for Figure 5)
omi_rankings_aw <- omi_resp |>
  mutate(pres_approval_dummy = case_when(
    pres_approval_dummy == 0 ~ "Putin critic",
    pres_approval_dummy == 1 ~ "Putin supporter"
  ),
  pres_approval_dummy = factor(pres_approval_dummy,
                               levels = c("Putin supporter",
                                          "Putin critic")),
  source_known_independent = case_when(
    source_known_independent == 1 ~ "Yes",
    source_known_independent == 0 ~ "No"
  ),
  source_known_independent = factor(source_known_independent,
                                          levels = c("Yes", "No"))) |>
  select(starts_with("ranking_"), "source_known_independent",
         "source_known_TV1", "source_known_Russia_1_24",
         "source_known_RIA", "source_known_RT", "source_known_RBC",
         "pres_approval_dummy", "female", "age_ordinal", "education_higher",
         "source_trusted_independent") |>
  select(-ends_with("Echo"), -ends_with("RBC")) |>
  mutate_at(vars(starts_with("ranking_")),
            list(~as.factor(.))) |>
  mutate_at(vars(starts_with("ranking_")),
            list(~relevel(., ref = "Don't know/hard to say"))) 

# dependent variables for the analysis (evaluations of state media outlets
# along various dimensions such as bias or accuracy)
ranking_list_aw <- omi_rankings_aw |> 
  select(starts_with("ranking")) |> names()

# run multinomial models (same as Figure 5, but conditional on knowledge
# of independent media)
rankings_regs_aw <- map(
  ranking_list_aw, 
  ~run_multinom(
    .x, reg_formula = " ~ pres_approval_dummy + source_known_independent + 
                        female + age_ordinal + education_higher +
                        pres_approval_dummy*source_known_independent",
    given_indep_knowledge = T))

# extract means (probabilities) of each category of evaluations
# conditional on presidential approval and knowledge of independent media
rankings_probs_aw <- map_df(rankings_regs_aw, "probs") |>
  filter(response != "Don.t.know.hard.to.say") |>
  mutate(eval_direction = case_when(
    response == "Mostly.yes" ~ "positive",
    response %in% c("Omit.important", "Often.gives.false.info", 
                    "Pro.government", "Anti.government", 
                    "Not.independent") ~ "negative",
  )) |>
  filter(eval_direction == "negative",
         pres_approval_dummy == "Putin supporter") |>
  mutate(criterion = case_when(
    dimension == "accuracy" & response == "Often.gives.false.info" ~ 
      "Often gives false info",
    dimension == "completeness" & response == "Omit.important" ~ 
      "Omits important facts",
    dimension == "bias" & response == "Pro.government" ~ 
      "Biased in favor of government",
    dimension == "independence" & response == "Not.independent" ~ 
      "Not independent"
  ),
  source = factor(source, levels = c("Channel One", "RIA", "RT", "Russia-24")),
  criterion = factor(criterion,
                     levels = c("Not independent",
                                "Biased in favor of government",
                                "Omits important facts",
                                "Often gives false info"))) |>
  # drop other dimensions (not needed for the figure)
  drop_na(criterion)

write_csv(rankings_probs_aw,
          "Estimates/Figure_D3_omi_perceptions_by_awareness_supporters.csv")

# Figure D3
ggplot(rankings_probs_aw, 
       aes(criterion, prob, 
           linetype = source_known_independent,
           shape = source_known_independent)) +
  geom_point(position = position_dodge(width = 0.6), size = 2) +
  geom_errorbar(aes(ymin = conf.low,
                    ymax = conf.high),
                width = 0,
                position = position_dodge(width = 0.6), size = 0.7) +
  labs(x = "", y = "Probability of agreeing with statements",
       title = "Perceptions of state media among Putin supporters",
       linetype = "Knows some independent media", 
       shape = "Knows some independent media") +
  ylim(c(0, 1)) +
  coord_flip() +
  facet_wrap(~source, ncol = 2) +
  theme_classic() +
  theme(legend.position = "bottom")
ggsave("Figures/Figure_D3.png")

# Additional APPENDIX Tables -------
# Table B1: Summary statistics ----------------

# summary statistics for the main study -----
resp_sum <- main_resp |>
  # limit the data set to respondents with non-missing presidential approval
  filter(session_id %in% main_stories$session_id,
         pres_approval_asked == 1) |>
  drop_na(pres_approval) |>
  select(pres_approval_dummy,
         source_independent_dummy,
         source_state_controlled_dummy, source_state_tv,
         female, education, age_group) |>
  mutate(age_group = case_when(
    age_group == "25-29" | age_group == "30-34" ~ "25-34", 
    age_group == "35-39" | age_group == "40-44" ~ "35-44", 
    age_group == "45-49" | age_group == "50-54" ~ "45-54", 
    age_group == "55-59" | age_group == "60-64" ~ "55-64", 
    TRUE ~ age_group 
  )) |>
  rename(`Approves of president (dummy)` = pres_approval_dummy,
         `Uses independent media` = source_independent_dummy,
         `Uses state-controlled media` = source_state_controlled_dummy,
         `Uses state TV` = source_state_tv,
         Female = female,
         `Higher education` = education,
         `Age` = age_group)

# create dummies for age groups
ag <- resp_sum |>
  to_dummy(`Age`, suffix = "label") |>
  rename_with(~str_replace(., "_", " "), everything())

# add dummies for age groups and summarize by subgroup
resp_sum <- resp_sum |>
  select(-`Age`) |>
  bind_cols(ag) |>
  summarise(across(everything(),
                   list(Mean = ~round(mean(., na.rm = T), 3),
                        N = ~sum(!is.na(.))))) |>
  pivot_longer(everything(), names_to = c("Variable", ".value"),
               names_sep = "_")

# summary statistics for the Levada study -----
levada_sum <- levada_resp |>
  # limit to respondents with non-missing presidential approval
  drop_na(putin_approval) |>
  mutate(`Uses state-controlled media` = NA_real_,
         `Uses independent media` = NA_real_) |>
  select(pres_approval_dummy,
         `Uses independent media`,
         `Uses state-controlled media`,
         info_federal_state_tv_dummy,
         female, education_higher, age) |>
  # create age groups analogous to the main study
  mutate(age_group = case_when(
           age <= 24 ~ "18-24",
           age %in% 25:34 ~ "25-34", 
           age %in% 35:44 ~ "35-44", 
           age %in% 45:54 ~ "45-54", 
           age %in% 55:64 ~ "55-64", 
           age > 64 ~ "65+" 
         )) |>
  select(-age) |>
  rename(`Approves of president (dummy)` = pres_approval_dummy,
         `Uses state TV` = info_federal_state_tv_dummy,
         Female = female,
         `Higher education` = education_higher,
         `Age` = age_group)

# create dummies for age groups
agl <- levada_sum |>
  to_dummy(`Age`, suffix = "label") |>
  rename_with(~str_replace(., "_", " "), everything())

# add dummies for age groups and summarize by subgroup
levada_sum <- levada_sum |>
  select(-`Age`) |>
  bind_cols(agl) |>
  summarise(across(everything(),
                   list(Mean = ~round(mean(., na.rm = T), 3),
                        N = ~sum(!is.na(.))))) |>
  pivot_longer(everything(), names_to = c("Variable", ".value"),
               names_sep = "_") |>
  mutate(Mean = if_else(Variable == "Uses state-controlled media" |
                          Variable == "Uses critical media",
                        NA_real_, Mean),
         N = if_else(Variable == "Uses state-controlled media" |
                       Variable == "Uses critical media",
                     NA_integer_, N))

# summary statistics for the OMI study ----
omi_sum <- omi_resp |>
  # limit to respondents with non-missing presidential approval
  drop_na(pres_approval) |>
  select(pres_approval_dummy,
         source_used_independent,
         source_used_state_controlled, source_used_state_tv,
         female, education_higher, age) |>
  rename(`Approves of president (dummy)` = pres_approval_dummy,
         `Uses critical media` = source_used_independent,
         `Uses state-controlled media` = source_used_state_controlled,
         `Uses state TV` = source_used_state_tv,
         Female = female,
         `Higher education` = education_higher,
         `Age` = age)

# create dummies for age groups
omi_ag <- omi_sum |>
  to_dummy(`Age`, suffix = "label") |>
  rename_with(~str_replace(., "_", " "), everything())

# add dummies for age groups and summarize by subgroup
omi_sum <- omi_sum |>
  select(-`Age`) |>
  bind_cols(omi_ag) |>
  summarise(across(everything(),
                   list(Mean = ~round(mean(., na.rm = T), 3),
                        N = ~sum(!is.na(.))))) |>
  pivot_longer(everything(), names_to = c("Variable", ".value"),
               names_sep = "_")

# combine summary stats for three studies into one table
# headers not fixed here; they are added in the paper directly
sumstats_resp <- bind_cols(resp_sum,
                           levada_sum |> select(-Variable),
                           omi_sum |> select(-Variable)) |>
  mutate(across(c(2, 4, 6), ~.x*100))

write_csv(sumstats_resp, "Estimates/Table_B1_sumstats_resp_3_surveys.csv")

# Table B2: News messages -------------------
# story characteristics and average level of belief in each news story

# baseline (sample-average) belief in each story (the proportion of respondents who
# found each story to be true)
main_beliefs_baseline <- main_stories |>
  group_by(story_code) |>
  summarise(Overall = mean(story_real)) |>
  ungroup() 

# belief in each story conditional on presidential approval 
# (the proportion of Putin supporters and critics who found each story to be true)
main_beliefs <- main_stories |>
  drop_na(pres_approval_dummy) |>
  mutate(pres_approval_dummy = if_else(pres_approval_dummy == 0,
                                       "Critic", "Supporter")) |>
  group_by(story_code, pres_approval_dummy) |>
  summarise(belief = mean(story_real)) |>
  ungroup() |> 
  pivot_wider(names_from = "pres_approval_dummy", values_from = "belief") |>
  # join with overall level of belief and story info/texts
  left_join(main_beliefs_baseline) |>
  left_join(main_stories_info |>
              select(story_code, political, story_direction, 
                     story_text_eng, fake)) |>
  # clean up
  arrange(story_code) |>
  mutate(across(c(Overall, Critic, Supporter),
                ~round(.x, 3)),
         Code = row_number(),
         `False?` = fake == 0,
         Political = if_else(political == 1, "Yes", "No")) |> 
  rename(Text = story_text_eng,
         Direction = story_direction) |>
  select(Code, Text, `False?`, Political, Direction, Overall, 
         Critic, Supporter)

# save estimates for the paper
write_csv(main_beliefs, "Estimates/Table_B2_story_summary_2020.csv")


# Table B3: Balance check -------------------
# chi-squared tests for equality of covariates across treatment groups
# calculated separately for each story, as treatments were assigned 
# independently for each story

# select columns from the story-level data set
main_stories_exp <- bind_cols(
  # sociodemographic covariates
  main_stories |>
    rename(Female = female,
           `Higher education` = education,
           `Age group` = age) |>
    select(story_label, story_source,
           Female,
           `Higher education`,
           `Age group`),
  # add columns for each subgroup by regime support
  main_stories |> 
    mutate(pres_approval_cat = case_when(
      pres_approval_cat == "Certainly approve" ~ "Strong supporter",
      pres_approval_cat == "Somewhat approve" ~ "Moderate supporter",
      pres_approval_cat == "Somewhat disapprove" ~ "Moderate critic",
      pres_approval_cat == "Certainly disapprove" ~ "Strong critic"
    )) |>
    to_dummy(pres_approval_cat, suffix = "label") |>
    rename_with(~str_remove(., "pres_approval_cat_")))

# extract covariate names for use in the function below
covariate_names <- names(main_stories_exp |> 
                           select(-story_label, -story_source))

# function to apply chi-squared tests to each column for each news story
compare_treatment_groups_chisq <- function(story) {
  
  # limit data to the selected story
  story_data <- main_stories_exp |>
    filter(story_label == story)
  
  # for each covariate, apply chi-squared test and extract p-value
  story_chisq <- map(covariate_names,
                     function(x) {
                       chisq_table <- story_data |>
                         drop_na({{ x }}) |>
                         select({{ x }}, story_source) |>
                         table()
                       story_df <- tibble(
                         covariate = x,
                         `News story` = story,
                         p_value = round(chisq.test(chisq_table)$p.value, 3)
                       )
                       return(story_df)
                     }) |>
    list_rbind()
  
  return(story_chisq)
}

# apply the function to all stories for all covariates of interest
assignment_by_covariates <- map(main_stories |> 
                                  filter(set == 1, 
                                         story_code > 103) |>
                                  distinct(story_label) |>
                                  pull(story_label),
                                compare_treatment_groups_chisq) |>
  list_rbind() |>
  # clean up
  arrange(`News story`) |>
  pivot_wider(names_from = covariate, values_from = p_value)|>
  relocate(`Age group`, .after = Female) |>
  relocate(`Strong supporter`, .after = `Higher education`) |>
  relocate(`Strong critic`, .after = `Strong supporter`) |>
  relocate(`Moderate critic`, .after = `Moderate supporter`)

# save calculations for the paper
write_rds(assignment_by_covariates, "Estimates/Table_B3_balance_check_main_study.rds")

